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ABSTRACT 

We have searched for continuous gravitational wave (CGW) signals produced by indi¬ 
vidually resolvable, circular supermassive black hole binaries (SMBHBs) in the latest EPTA 
dataset, which consists of ultra-precise timing data on 41 millisecond pulsars. We develop fre- 
quentist and Bayesian detection algorithms to search both for monochromatic and frequency- 
evolving systems. None of the adopted algorithms show evidence for the presence of such a 
CGW signal, indicating that the data are best described by pulsar and radiometer noise only. 
Depending on the adopted detection algorithm, the 95% upper limit on the sky-averaged strain 
amplitude lies in the range 6 x 10“ < A < 1.5 X 10 at 5nHz < / < 7nHz. This limit 

varies by a factor of five, depending on the assumed source position, and the most constraining 
limit is achieved towards the positions of the most sensitive pulsars in the timing array. The 
most robust upper limit - obtained via a full Bayesian analysis searching simultaneously over 
the signal and pulsar noise on the subset of ours six best pulsars - is A « 10“^"^. These limits, 
the most stringent to date at / < lOnHz, exclude the presence of sub-centiparsec binaries 
with chirp mass Ale > 10 ®Mq out to a distance of about 25Mpc, and with Me > IO^^Mq 
out to a distance of about IGpc {z ~ 0.2). We show that state-of-the-art SMBHB population 
models predict < 1% probability of detecting a CGW with the current EPTA dataset, consis¬ 
tent with the reported non-detection. We stress, however, that PTA limits on individual CGW 
have improved by almost an order of magnitude in the last five years. The continuing advances 
in pulsar timing data acquisition and analysis techniques will allow for strong astrophysical 
© 2010 RAS constraints on the population of nearby SMBHBs in the coming years. 
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1 INTRODUCTION 


The direct detection of gravitational waves (GWs) is one of the pri¬ 
mary goals of contemporary observational astrophysics. The access 
to GW information alongside well established electromagnetic ob¬ 
servations will be a milestone in our investigation of the Universe, 
opening the era of multimessenger astronomy. 

Precision timing of an array of millisecond pulsars (i.e. a pul¬ 
sar timing array, PTA) provides a unique opportunity to get the 
very first low-frequency (nHz) GW detection. PTAs exploit the ef¬ 
fect of GWs on the propagation of radio signals from ultra-stable 
millisecond pulsars (MSPs) to the Earth (e.g. |Sazhin|[T978[ |De-| 
|tweiler|[l979^ , producing a characteristic fingerprint in the times 
of arrival (TOAs) of radio pulses. In the timing analysis, TOAs 
are fitted to a physical model accounting for all the known pro¬ 
cesses affecting the generation, propagation and detection of the 
radio pulses. The timing residuals are the difference between the 
observed TOAs and the TOAs predicted by the best-fit model, 
and they carry information about unaccounted noise and poten¬ 
tially unmodelled physical effects, such as GWs, in the datas- 
tream (e.g. |Hellings & Downs] |1983| |Jenet et al.|[ 


s|[T^JT 

ropean Pulsar Timing Array (EPTA, [Kramer & Champion|2013| >, 
the Parkes Pulsar Timing Array (PPTA, |Hobbs||2013^ and file 
North American Nanohertz Observatory for Gravitational Waves 
(NANOGrav, |McLa ughlin|20 1 3| l, joining together in the Interna¬ 
tional Pulsar Timing Array (IPTA, |Hobbs et al.|2010[|Manchesf^ 
|& IPTA|2013| l, are constantly improving their sensitivity in the fre¬ 
quency range of ~ 10“® — 10“® Hz. 

The primary GW source in the nHz window is a large popula¬ 
tion of adiabatically inspiralling supermassive black hole binaries 
(SMBHBs), formed following the frequent galaxy mergers occur¬ 
ring in the Universe jBegelman et al.|1980^. Signals from a cosmic 
string network (see, e.g. |Vilenkin|1981[|Vilenkin & Shellard|1994) 
or from other physical processes occurring in the early Universe 
(see, e.g. |Grishchuk|2005) are also possible, but we will concentrate 
on SMBHBs in this paper. Consisting of a superposition of sev¬ 
eral thousands of sources randomly distributed over the sky ( |Sesana| 
|et al.|200^ , the signal has classically been described as a stochas¬ 
tic GW background (GWB, [Rajagopal & Romani||1995[ |Jaffe &| 


[2005). 


The Eu- 


|Backer||2003[ [Wyithe & Loeb||2003[|Sesana et al.|2004| . Conse¬ 

quently, in the last decade several detection techniques have been 
developed in this direction (e.g. Anholm et al. 2009 van Haasterenj 


jet al.|2009[|Lentati et al.|2013||Chamberlin et al.|2014) and applied 

to the EPTA, PPTA, and NANOGrav datasets to get limits on the 
amplitude of a putative isotropic GWB ( Jene£^£^alJ2006 JTaid^ 


et al.|20lT]|van Haasteren et al.|201 l||Demorest et al.|2013[|Sh' 


non et al.|2013[|Lentati et al.|2015| l. 

However, ISesana et al.|P009l > (see also |Ravi et al.|201^ first 
showed that the signal is dominated by a handful of sources, some 
of which might be individually resolvable. The typical evolution 
timescale of those SMBHBs is thousands-to-millions of years, far 
exceeding the observational baseline of PTA experiments (about 
two decades); therefore, their signals can be modeled as non¬ 
evolving continuous GWs (CGWs, [Sesana & Vecchio|2010) . Re¬ 
solvable sources are particularly appealing because, if detected and 
localized on the sky, they can also be followed up electromagnet- 
ically, thus providing a multimessenger view ( Sesana^etjil J2^2) 


[Tanaka et al.|201^|Burke-Spolaor|2013[|^sado & Sesana|2014 1 . 

This prospect triggered a burst of activity in the development 
of search and parameter estimation algorithms for CGWs from cir¬ 
cular SMBHBs ([S 


& Vecchio|20T0l|Lee et al.|2011 


|& Sesana||20T2l jEllis et al.||20T2l jPetiteau et ar]|20T3a 


Babak 


Taylor 


jet al.|20r4) , and more recently led to the development of the first 
pipelines for eccentric binaries (Tnylor etay2015[ |Zhu et al.|2015) . 
The pioneering work of jYardley et al. ^2010^ was the first to pro¬ 
duce sensitivity curves and set upper limits using a power spectral 
summation method. More recently, [Arzoumanian et ar] ( |2014) ap¬ 
plied the frequentist and Bayesian methods for evolving and non¬ 
evolving signals described in Ellis et aL[ )2012) to the NANOGrav 
5-year dataset ( [Demorest et al.||2013) , whereas [Zhu et ak[ j2014) 
applied a frequentist method to the PPTA data release (DRl) pre¬ 
sented in [Manchester et ar|[2013) . Those limits are usually cast 
in terms of the intrinsic strain amplitude of the wave, ho or its 
inclination-averaged version (which is a factor of 1.26 larger) as 
a function of frequency, both averaged over the entire sky or as 
a function of sky location. The best sky-averaged 95% c onfidence 
upper limit on ho quoted to date is 1.7 x at lOnHz i Zhu et al. 
[ 2 ^ . 

Here we investigate the presence of non-evolving continuous 
waves from circular binaries in the latest EPTA data release ( |Desvi-| 
). We perform a comprehensive study applying 


both frequentist dBabak & Sesana 

2012 [Ellis et al. [2012[ jPetiteau 

et al.[|2013ai and Bayesian (Ellis 

20 13[ [Taylor et al ||2014[ Las- 


sus et al.|2015) methods, and searching for both evolving and non¬ 


evolving GW signals. The paper is organized as follows. In Sec- 
tion[^we introduce the EPTA dataset and the adopted gravitational 
waveform model. Section [^ is devoted to the description of the 
techniques developed to analyze the data, divided into frequentist 
and Bayesian methods. Our main results (upper limit, sensitivity 
curves, sky maps) are presented in Section [^ and their astrophysi- 
cal interpretation is discussed in Section[^ Finally, we summarize 
our study in Section [^ Throughout the paper we use geometrical 
units G = c = 1. 

This research is the result of the common effort to directly 
detect gravitational waves using pulsar timing, known as the Euro¬ 
pean Pulsar Timing Array (Kramer & Champion 2013^ 


2 EPTA DATASET AND GRAVITATIONAL WAVE 
MODEL 

2.1 The EPTA Dataset 

In this paper we make use of the full EPTA data release described 
in |Desvignes et alT[ ( |2015) , which consists of 42 MSPs monitored 
for timespans ranging from 7 to 24 years. However, we exclude 
PSR 11939-1-2134 from our analysis because it shows a large, un¬ 
modelled red-noise component in its timing residuals. The remain¬ 
ing 41 MSPs show well-behaved rms residuals between 130ns and 
35/is. For each of these pulsars, a full timing analysis has been 
performed using a time-domain Bayesian method based on Multi- 
Nest ( [Feroz et al.|2009) , which simultaneously includes the white 
noise modifiers EFAC and EQUAE0 for each observing system, 
as well as intrinsic red noise and (observational) frequency depen¬ 
dent dispersion measure (DM) variations. Variations in the DM are 
due to a changing line-of-sight through the interstellar medium to¬ 
wards the pulsar. Hereinafter, we refer to this timing analysis as 
single pulsar analysis (SPA). As a sanity check, parallel analyses 


1 www.epta.eu.org/ 

^ EFAC is used to account for possible mis-calibration of the radiometer 
noise and it acts as a multiplier for all the TOA error bars. EQUAD repre¬ 
sents some additional (unaccounted) source of time independent noise and 
it is added in quadratures to the TOA error bars. 
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have also been done in the frequency domain using the TempoNest 
plugin I Lentelietalj 20T4l l for the Tempo2 pulsar timing pack¬ 
age ob^r’eT^ |2006^ , and in time domain using search method 
combining a genetic algorithm ( Petiteau et al.|201^ with MCM- 
CHammer ( [Foreman-Mackey et al.|2013| ). The results of the three 
methodologies are consistent. 

For the searches performed in this paper, we use the results 
of the SPA produced by MultiNest. Those consist of posterior 
probability distributions for the relevant noise parameters (EFAC, 
EQUAD, DM and intrinsic red noise), together with their max¬ 
imum likelihood (ML) values. Extensive noise analyses on the 
same dataset are fully detailed in |Janssen et al.| (| 2015^ ; [Uaballero| 
|et al.| ^2015^ , and the posterior distributions of the noise parame¬ 
ters of the three most sensitive pulsars in our array (J1909—3744, 
J1713+0747, J1744—1134) are also given in Figure 3 of the com¬ 
panion paper |Lentati et al.| P015| l. In most cases we fix the 
noise parameters to their ML value, but we also perform separate 
searches sampling from the posterior distribution or searching si¬ 
multaneously over the signal and the noise parameters, as described 
in Section[3and summarized in Table[T] 


2.2 Gravitational wave model 

In this subsection we introduce the mathematical description of the 
GW signal from a binary in circular orbit and the associated PTA 
response. The main purpose here is to introduce the notation; we 
refer the reader to |Sesana & Vecchio| ( |2010t for a full derivation 
of the relevant equations. The timing residuals of radio pulses due 
to the propagation of the electromagnetic waves in the field of an 
intervening GW can be written as 


where 

fa(f) = — {(1 + cos^ i)F^ [sin(u)f -|- <l>o) - sin<Fo] + 

OJ 

2 cos iFa [cos(ajf -I- $ 0 ) — cos (Do]} , 

fa(f) = — {(1 + COS^ b)F^ [sin(u;af -|- -Fa -b 4)0)- 

sin(<Fa -b $ 0 )] + 2 cos lF^ [cos(cUaf -b 'Fa -b $ 0 )- 

cos(‘Fa -b 43 o)]} . 

( 6 ) 

Equation contains all the relevant features of the signal, and 
deserves some detailed explanation. Here b is the inclination of 
the SMBHB orbital plane with respect to the line-of-sight, and A 
(sometimes referred to as ho in the PTA literature) is the amplitude 
of the GW strain given by 

(7) 

Throughout the paper we consider a SMBHB with redshifted chirp 
mass A4c = where M = (mi-bm 2 ) and r] = mim 2 /M^ 

are the total mass and symmetric mass ratio of the binary system; 
Dl is the luminosity distance to the source, and / = u;/(27r) is 
the observed GW frequency, which is twice the SMBHB observed 
orbital frequenc)]^ We specify A and Aa in equation because 
the GW frequency might be different in the pulsar and Earth terms, 
implying different amplitudes. In fact, in the quadrupole approxi¬ 
mation, the evolution of the binary orbital frequency oJorh = a;/2 
and GW phase can be written as 

CCorb(f) = tUorb , (8) 


where 



StVg 

Va 


1 faPj 

2 1 -bPa-G 


A/ii 


(I) 


( 2 ) 


Here Vg is the frequency of the TOAs (i.e., the spin frequency of 
the pulsar) and Svg its deviation. The index a labels a particular 
pulsar (a = 1, ...,N where A = 41 is the number of the pulsars 
in our array) and indicates that a given quantity explicitly depends 
on the parameters of the individual pulsar, p denotes the position of 
the pulsar on the sky, and Cl is the direction of the GW propagation. 
The last factor in equation depends on the strain of the GW at 
the location of the pulsar hij{tg) and on Earth hij{t): 


Ahij = hij(tl) - hij{t)- (3) 


^ • (9) 

In equation oj and ujg = utit — Tg) are the GW frequencies of 
the Earth term and pulsar term, respectively. Over the typical dura¬ 
tion of a PTA experiment (decades), these two frequencies can be 
approximated as constants (see, e.g. ISesana & Vecchio|2010) , and 
we drop the time dependence accordingly. However, the delay Tg 
between the pulsar and the Earth term is of the order of the pulsar- 
Earth light travel time, and can be thousands of years. This is com¬ 
parable with the evolutionary timescale of the SMBHB’s orbital 
frequency ^Sesana & Vecchio|2010| >. In particular u)g = 2u;orb,a 
is obtained by setting t = —Ta in the right-hand side of equa¬ 
tion and can be generally different from one pulsar to the other 
and from uj. Combining equations and l|^, one can show that 

"Fa « -i^Tg. 

If T is the timing baseline of a given PTA’s set of observations, 
then its nominal Fourier frequency resolution bin is given by A/ « 
1/T. We therefore have two possibilities: 


The pulsar time is related to the Earth time t as: 


• If {ujg — ljj)/{2'k) > A/ for the majority of the MSPs in the 


Fg = t - Lg{l + Cl.pg) =t - Tg, (4) 

where Lg is the distance to the pulsar. We consider a non-spinning 
binary system in quasi-circular orbit. To leading order, the response 
of a particular pulsar to a passing GW (that is, the induced timing 
residual) is given by: 

ra{t) = r^{t) (5) 


^ The relation between redshifted chirp mass and rest-frame chirp mass is 
j\4c = (1 + z)J\Ac,rf’ likewise, the relation between observed and rest- 
frame frequency is /^f = (1 + 2 :)/. The GW community works with red¬ 
shifted quantities because this is what is directly measured in the observa¬ 
tions, and because this way, the 1 + 2 : factors cancel out in the equations, 
making the math cleaner. We note, moreover, that current PTAs might be 
able to resolve SMBHBs out to 2 : ~ 0.1, implying only a small difference 
between rest-frame and redshifted quantities. 
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array, then the pulsar and the Earth terms fall in different frequency 
bins; all the Earth terms can be added-up coherently and the pulsar 
terms can be considered either as separate components of signal or 
as an extra incoherent source of noise. 

• Conversely, if {uja — uj)/{2-k) < Af for the majority of the 
MSPs, then the pulsar terms add up to the respective Earth terms, 
destroying their phase coherency. 

This distinction has an impact on the detection methodology that 
should be adopted, as we will see in the next Section. 

The antenna response functions of each pulsar to the 

GW signal depend on the mutual pulsar-source position in the sky 
and are given by: 



1 {r.uf - 

(10) 


2 1 + p^-.Ci 

= 

{p°’.u){p°’.v) 

(11) 


1 -1- p“.D 


where 



D = 

— {sin 9s cos ps, sin 9s sin ps , cos 0s }, 

(12) 

u — 

h cos Ip — m sin Ip, v = mcosp + nsinip. 

(13) 

h = 

(cos 9s cos ps , cos 0s sin ps , — sin 0s }, 

(14) 

771 = 

(sin ps,- cos 0s, 0}. 

(15) 


Here 6s,(i)s define the source sky location (respectively latitude 
and longitude), and ip is the GW polarization angle. 


2.3 Likelihood function and noise model 


The core aspect of all searches performed in this study (both fre- 
quentist and Bayesian) is the evaluation of the likelihood that some 
signal described by equation is present in the time series of the 
pulsar TOAs. We use the expression for the likelihood marginalised 
over the timing parameters as described in detail in | van Haasteren| 
|& LevinldM^ : 


exp 


^(27r)—def(G^CG) 

- ^'^G{G^CG)-^G^{5t - r) 


( 16 ) 


Here n is the length of the vector 5t = UXa obtained by concate¬ 
nating the individual pulsar TOA series Xa, m is the total num¬ 
ber of parameters describing the timing model, and the matrix G 
is related to the design matrix (see [van Haasteren & Levin|2013| 
for details). The variance-covariance matrix G, in its more gen¬ 
eral version, contains contributions from the GWB and from the 
white and (in general) red noise: G = Cgw -f G^n -f Gm- We 
I Iva n Haasteren et^p009] l and to our compan- 

noise 


refer the reader to ' _ 

paper [Lentati et al. ( 2015^ for exact expressions of the : 


ion 


variance-covariance matrix. In our analysis we use both the time 
d vai^aasteren & Levin|20T^ and frequency domain ( |Lentati et al.| 
|2013 1 representation of this matrix. Both approaches give qualita¬ 
tively and quantitatively consistent results (as we checked during 
our analysis). Therefore, we will not specify which representation 
was used for each of the employed methods. Moreover, we have 
excluded Cgw from our analysis assuming the hypothesis that the 
data consists of noise and a single CGW source. 

The model parameters in equation jl 6 | l are split in two groups 

(i) parameters describing the CGW signal (A : r = r(A)), and 

(ii) parameters describing the noise model 9. The waveform of a 


non-spinning circular binary given by equation 0 is generally de¬ 
scribed by 7 + 2N parameters. The Earth term (a single sinusoidal 
GW) is fully described by 7 parameters: (.4, 9s, (ps, T*, t, w, $ 0 ), 
and each pulsar term adds two additional parameters: frequency 
and phase [uJa, ‘f’a). As mentioned above, the pulsar term might 
fall at the same frequency as the Earth term, in which case we have 
only one extra parameter per pulsar (since uja = io). In principle, 
even for uja 7 ^ a;, Ua and <l?a are uniquely connected through the 
pulsar distance La as shown by equations and 1 ^. However this 
implicitly assumes that we have an exact model for the evolution of 
the binary, which in this case is perfectly circular, non-spinning, 
and GW-driven. Even tiny deviations from these assumptions (i.e., 
small residual eccentricity, partial coupling with the environment, 
spins), very likely to occur in nature, will invalidate the uia — 
connection, and in the most general case, the two parameters must 
be considered separately. We will specify the details of the adopted 
waveform model for each individual method in the next section. 

Some of the noise parameters (like pulsar-intrinsic red noise 
or clock and ephemeris errors) are correlated with the GW param¬ 
eters and one should in principle fit for noise and GW parameters 
simultaneously. However, such a multidimensional search is com¬ 
putationally very expensive, and in most of the searches detailed 
below, we use noise characteristics derived from the SPA intro- 
and extensively described in|Janssen et al.| 


duced in Section 


2.1 


( |2015^ ; |Caballero et al.H2015^ . We characterise the noise by con¬ 
sidering the data from each pulsar separately (as given by the 
SPA), and to exclude any potential bias we also considered a model 
“noise - 1 - monochromatic signal”. The results of the latter are usu¬ 
ally consistent with the “noise only” model, except for one pul¬ 
sar , J1713-I-0747, where we have found some correlation between 
the parameters describing the red noise and the extra monochro¬ 
matic signal. Each SPA returns posterior distributions for the array 
of noise parameters 9a’. slope and amplitude of the red noise, slope 
and amplitude of DM variations (both red noise and DM varia¬ 
tion were modeled as single power laws) and EFAC and EQUAD 
for each backend. We have used this information in two ways: (i) 
use the ML estimator for all the parameters ( 9ml = ) 

and assume that the noise is represented by that model, (ii) draw 
the parameters describing the noise from the posterior distributions 
obtained in the SPA. The first choice (fixed noise at 0 ml) is compu¬ 
tationally cheaper as we need to compute the noise variance matrix, 
C{9), only once, while in the second choice we need to recompute 
it for each draw of the noise parameters. [Arzoumanian et al.lpOl^ 
found that fixing the noise to the ML values will bias the results 
of the search and will yield a better upper limit compared to the 
full search including noise parameters. The effect is, however, not 
dramatic, as they found their upper limit is less than a factor « 1.5 
worse in the latter case, over the full frequency range. To test the ro¬ 
bustness of our results we also ran a full analysis searching simulta¬ 
neously on the GW and noise parameters on a restricted dataset of 6 
pulsars (see Section[T3]for the definition of this restricted dataset). 


3 METHODS 

As described in the previous section, the two parts of the GW sig¬ 
nal (Earth term and pulsar term) might or might not fall at the 
same frequency, which has implications for the form of the like¬ 
lihood function given in equation On top of that, both fre- 
quentist and Bayesian methods can be used to analyze the data. 
We therefore identify four separate classes of analysis: frequen- 
tist non-evolving, frequentist evolving, Bayesian non-evolving, and 
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Bayesian evolving. In each individual case, the signal (and some¬ 
times also the noise) can be treated differently, and the likelihood 
might undergo peculiar manipulations (maximization, marginaliza¬ 
tion, etc.); moreover, we sometimes explore two variations of the 
same method. The result is a variety of complementary searches, 
which we now describe in detail and which are comparatively sum¬ 
marized in Table [T] 


With the assumption of Gaussian noise, 2Tp follows a cen¬ 
tral distribution, po{J-p) - non-central, p\{Tp,p), if a signal 
with optimal S/N = p is present - with n = 2N degrees of free¬ 
dom (which follows from 2N maximisations of the log likelihood), 
given by: 

-pn/2-\ 


3.1 Frequentist methods 

The basis for the frequentist analysis is to postulate that we have 
a deterministic signal in the data which is either corrupted (in case 
of detection) or dominated by noise (no detection). We then define 
appropriate statistical distributions (or simply “statistics”) based on 
the likelihood function, those describe the data in absence and pres¬ 
ence of a signal. Those statistics must be chosen in such a way that 
the detection rate is maximised for a fixed value of the false alarm 
probability (TAP), which is also known as the Neyman-Pearson cri¬ 
terion. The aim is to check the null hypothesis (i.e., whether the 
data are described by noise only), and in case of “no detection” 
set an upper limit on the GW amplitude. In building our statistics 
we always assume that the noise is Gaussian and, unless otherwise 
stated, characterised by ML parameters estimated during the SPA. 
We then fix the FAP at 1% to set the detection threshold. In case 
of no detection, frequency dependent upper limits are obtained by 
splitting the frequency range in small bins, performing a large num¬ 
ber of signal injections with varying amplitude in each bin, and 
computing the associated detection statistics. In the next two sub¬ 
sections we describe two particular implementations of this proce¬ 
dure, known as tFp and Te. statistics. 


3.1.1 Tp-statistic 

In the case of non-evolving sources (i.e. uja = w), |Ellis et ^ 
( |2012^ showed that, for each pulsar a, one can write equation 
factorizing out the cu dependence 






Explicit expressions for and kq) can be found in Ellis et al. 


1 2012|. We merely stress here that are independent on the con¬ 
sidered pulsar a. The log of expression l |16^ can then be maximized 
over the constant 2N coefficients fo(i,a), b{ 2 ,a) assuming that they 
are independent, resulting in what is commonly known as the IFp- 
statistic. This latter assumption is not actually true if the number of 
pulsars is larger than 6. We make use of the full 41 pulsar dataset 
in the J^p-statistic evaluation. However, 6 pulsars dominate (i.e., 
give 90% of) the S/N, as discussed in Section [X3| moreover, we 
can simply postulate this form of statistic. Assuming independence 
in the maximisation process will somewhat degrade the detection 
power of the J^p-statistic; nonetheless, this is a very simple detec¬ 
tion statistic which depends only on one parameter: the frequency 
of the GW. The J^p-statistic is in essence an “excess-power” statis¬ 
tic, in which we basically search for extra power -compared to 
the expectations of the statistic describing the null hypothesis- at 
a given frequency in each pulsar’s data (summing the weighted 
square of the Eourier amplitudes in each pulsar). By mixing the 
GW phases <I?o and <!>„ at the Earth and at the pulsar in the coeffi¬ 
cients by construction IFp assumes that there is no coherence 
between GW signals across each individual pulsar’s data. 




^n/2—1 


(19) 


where I„/ 2 -i{x) is the modified Bessel function of the first kind 
of order n/2 — 1. We divide our dataset in bins A/ = IjT, where 
T is the longest pulsar observation time, and we evaluate IFp in¬ 
dependently in each bin. At any given /, we consider only pulsars 
with observation time T > 1/ f when evaluating IFp. Since the im¬ 
plementation of the -statistic is computationally cheap, we run 
two flavors of it: (i) one in which we fix the noise to the ML values 
(J-p-ML hereinafter), (ii) one in which we take into account the un¬ 
certainty in the noise parameters by sampling from their posterior 
distribution derived from the SPA (simply IFp hereinafter). 


3.1.2 IFe-statislic 


If the source is evolving, then oJa ^ io and one can choose to con¬ 
sider as “signal” only the Earth portion of equation (j^, and treat 
the pulsar term as a source of noise. In this case, we take only the 
Earth term of r in the likelihood expression l |16[ l. |Babak & Sesana| 
l |2012^ showed that it is convenient to rewrite the antenna pattern 
expressions (|TTJ separating the terms containing the polarization 
angle ip: 


F+ = F: cos{2iP) + F: sm{2pj) 

FPP = -F;cos(2t/>)+F“sin(2t/.), 

thus rearranging equation in the form 


( 20 ) 

( 21 ) 


Xait) = ( 22 ) 




Babak 


Explicit expressions for ayj, F^ can be found in 

|& Sesan^ ( |2012^ , but note that, contrary to the Fp case, the coetti 
cients a(j) do not depend on the considered pulsar a. As done for 
the Fp case, we can now maximize over the four constants a(i), 
constructing a statistic - the J^e-statistic - which is a function of 
three parameters; namely the GW frequency, /, and source sky lo¬ 
cation, 9s,<ps- If the noise is Gaussian, Fe follows a x^ distri¬ 
bution, ppi{FA - non-central, p\{F^,p), if a signal with optimal 
S/N= p is present - with four degrees of freedom (which follows 
from the four maximisations of the log likelihood), given by: 


Po{Fe) = Fee 


(23) 


px{Fe,p) = (24) 

Where /i(a:) is the modified Bessel function of the first kind of 
order 1. 

Note that the J^e-statistic adds the data from the pulsars co¬ 
herently (taking the phase information into account), and being a 
function of the source position, it allows direct sky localization. 
However, if the signal is non-evolving, its efficiency drops signifi¬ 
cantly and the estimation of the source sky location can be severely 
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biased. The evaluation of is significantly more computationally 
expensive than as we need to take into account cross-pulsar 
correlations. We therefore continue to use the full 41-pulsar dataset, 
but we fix the noise parameters at the ML values found in the SPA, 
without attempting any sampling (contrary to the case). 

For searching over the 3-D parameter space (/, Q, 0), we 
use the multi-modal genetic algorithm described in |Petiteau et al.| 
( |2013b^ . A detailed description of the method can be found there; 
here we give only a brief overview. In the first step we run a ge¬ 
netic algorithm (with 64 organisms) over 1000 generations tuned 
for an efficient exploration of the whole parameter space. Then we 
identify the best spot in the 3-D space and seed there a variation 
of the Markov chain Monte Carlo "MCMC Hammer" ( |Foremam| 
[Mackey et al.|201^ , which serves as a sampler and returns the ef¬ 
fective size of the “mode” and the correlations among parameters. 
Here “mode” stands for a local maximum of the likelihood. We then 
again run the genetic algorithm, exploring the remaining parame¬ 
ter space (i.e. excluding the mode we already found) and searching 
for other potential local maxima in the likelihood. At each distinct 
maximum that is found, we run "MCMC hammer", and we iterate 
this procedure 5 times. The end result is a set of independent local 
maxima in likelihood, among which we choose the largest (high¬ 
est in value). We repeat the entire procedure three times to verify 
convergence. This algorithm has proven to give very robust results, 
even in the case of pathological likelihood surfaces with multiple 
maxima. 


3.2 Bayesian analysis 


In the Bayesian approach, the parameters describing the model are 
treated as random variables. The principles of Bayesian statistics 
provide a robust framework to obtain probability distributions of 
those parameters for a given set of observations, while also folding 
in our prior knowledge of them. 

Bayes’ theorem states that the posterior probability density 
function (PDF), p(/2|d, H), of the parameters /2 within a hypothe¬ 
sised model T-i and given data d is 


p{fl\d,'H) = 


£(/2,'H|d)7r(/2|'H) 

p(djW) 


(25) 


where C{fi, 'H\d) is the likelihood of the parameters given the data 
assuming the model T-i with parameters fl. The prior PDF of the 
model parameters, 7r(/2|'H), incorporates our preconceptions and 
prior experience of the parameter space. The Bayesian evidence, 
p{d\'H) = Z, is the probability of the observed data given the 
model %, 




Z = j C{fT)Tt{fl)d^p,. 


(26) 


For posterior inference within a model, Z plays the role of a nor¬ 
malisation constant and can be ignored. However, if we want to 
perform model selection then the evidence value becomes a key 
component in the computation of the posterior odds ratio: 


p{'H2\d) _ p{d\'H2)pi'H2) _ Z 2 X p{H2) 


(27) 


p{Hi\d) p{d\'Hi)p{Hx) ^ixp(-Hi)' 

Here Tfi, 'H 2 are the two models under comparison, Z^jZ\ is the 
Bayes factor, andp('H2)/p(7fi) is the prior probability ratio for the 
two competing models, which is often set to be one; the posterior 
odds ratio is then just the Bayes factor. 

When we specialize the formalism above to PTA data, the like¬ 
lihood L is given by equation 1 16 s the data d is the concatenation of 


the TOA series and the model’s parameters fl are identified with 
9, A. The presence of a signal in the data is assessed through the 
odds ratio, equation l |27^ , where model Hi corresponds to data de¬ 
scribed by noise only and model H 2 corresponds to data described 
by noise plus a CGW. We use non-informative and conservative 
priors in our analysis: they are always uniform in all parameters. It 
is especially important to emphasise that we have used a broad uni¬ 
form prior in the signal’s amplitude A, which will provide robust 
and conservative upper limits on the strain. Unless otherwise stated, 
we fix the stochastic noise parameters to the ML values found in 
the SPA, and we employ either Multi Nest (jFeroz e t al.||2009^ or 
a parallel tempering adaptive MCMC ( [Lassus et al.||2015^ to re¬ 
turn samples from the posterior distributions, and thus reconstruct 
the posterior PDFs. Both techniques also permit a recovery of the 
aforementioned model-comparison statistic known as the Bayesian 
evidence. We describe the different searches in detail below. 


3.2.1 Phase-marginalised Bayesian analysis 

For the non-evolving CGW signal searches, one should sample a 
7 + N multidimensional posterior, corresponding to the parameters 
For a 41 -pulsar dataset, that amounts 
to sampling a 48-dimensional parameter space. We mitigate the 
computational cost implied by such high-dimensionality by numer¬ 
ically marginalising (integrating) our posterior distribution over all 
of these nuisance parameters <I>a, thereby collapsing the search to 
more manageable dimensions (i.e. to seven parameters only,|Taylor| 
[et al.|2014[ (. By doing so, we not only rapidly recover the poste¬ 
rior PDFs, but also achieve quick and accurate Bayesian evidence 
values. This method is close in spirit to the J^p-statistic, we call it 
phase-marginalised Bayesian analysis (labeled as Bayes-EP-NoEv, 
for Bayes, Earth-l-pulsar, non-evolving source). The sampling of the 
posterior is performed by MultiNest. 

The high performance of the MultiNest sampler allowed us to 
also run a full search, including noise parameters 6, on a restricted 
dataset composed of the 6 best pulsars -contributing to 90% of the 
total S/N, see Section [T3] and figure[T]- in our array (labeled Bayes- 
EP-NoEv-noise). We use non-informative priors also for the pulsar 
noise parameters: uniform in the range [1,7] for the slopes of both 
red noise and DM; uniform in the range [—20, —10] for the log¬ 
arithm (base 10) of their amplitudes; uniform in the range [0,10] 
for the global EFAC (see [Lentati et al.|[2015] for a definition of 
global EEAC). The posterior spans now a 7 -F 1>N — 37 dimen¬ 
sion parameter space. We designed this search as a benchmark for 
our different fixed-noise analysis; to speed up the sampling, we re¬ 
stricted the frequency prior range to [5,15]nHz, which turns out to 
be the sweet-spot of our array sensitivity. 


3.2.2 Full Bayesian analysis 

We also employed the Bayesian formalism to construct a search 
for evolving signals. In this case we can either (i) use the full sig¬ 
nal of equation or (ii) restrict ourselves to the Earth term only. 
This latter case (ii) is similar in spirit to the frequentist J7=-statistic, 
with the difference that we now search over the full 7-dimensional 
source parameter space with the whole 41-pulsar array; we label 
this analysis Bayes-E. The full-signal analysis (i) is very compu¬ 
tationally expensive because we do not assume any model for the 
source evolution, therefore adding two extra parameters (uia, $a) 
for each pulsar. Note that our parametrisation of the full signal 
does not rely on any knowledge about the distance to the pulsar. 
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Search ID 

Noise treatment 

pulsars (A^) 

parameters 

Signal model 

Likelihood 

Fp-ML 

Fixed ML 

41 

1 

E+P NoEv 

Maximization over 2N constant amplitudes 

Fp 

Sampling posterior 

41 

1 

E+P NoEv 

Maximization over 27V constant amplitudes 

Fe 

Fixed ML 

41 

3 

E 

Maximization over 4 constant amplitudes 

Bayes_E 

Fixed ML 

41 

7 

E 

Full 

Bayes_EP 

Fixed ML 

6 

7 + 2x6 

E+PEv 

Full 

Bayes_EP_NoEv 

Fixed ML 

41 

7 

E+P NoEv 

Marginalization over pulsar phases 

Bayes_EP_NoEv_noise 

Seai'ched over 

6 

7 + 5x6 

E+P NoEv 

Marginalization over pulsar phases 


Table 1. Summary of the searches performed in this study. Column 1: name of the search; column 2: treatment of the noise in the search; column 3; number of 
pulsars considered in the dataset (Ai); column 4; dimensionality of the parameter space to search over; column 5; adopted signal model; column 6; notes about 
the treatment of the likelihood function. The different signal models are; Earth term only (E), Earth plus pulsar term at the the same frequency (E+P NoEv), 
Earth plus pulsar term at different frequencies (E+P Ev). 



Figure 1. Cumulative normalised (S/N)^. We rank pulsars according to their relative contribution to the total S/N, and we plot the quantity E,i<ca P°-lP- 


We name this search Bayes-EP-Ev. We use non-informative pri¬ 
ors also on those parameters, imposing the only constraint that the 
pulsar-term frequencies cannot be higher than the Earth-term one 
(uJa SS cj). Since we need to cover a 7 + 2N parameter space, we 
limit this search to the best 6 pulsars, for a total of 19 parameters. 
In both searches, the multidimensional posteriors are stochastically 
sampled with a parallel tempering adaptive MCMC (|Lassus et al.| 
|20T5l l. 


3.3 Restricted dataset: ranking pulsar residuals 

As mentioned above, some of the searches are extremely compu¬ 
tationally expensive, involving sampling of highly dimensional pa¬ 
rameter spaces. A way to mitigate the computational cost is to run 
the algorithms on a “restricted dataset”, which includes only the 
best pulsars for these purposes. We therefore need a metric to rank 
the quality of each pulsar. We choose our metric to be the rela¬ 


tive S/N contribution of each pulsar to a putative detectable source. 
We conduct extensive Monte Carlo simulations, in which we inject 
CGW signals with an astrophysically motivated distribution of pa¬ 
rameters A in the EPTA dataset. For each signal we compute the 
total S/N and the relative contribution of each individual pulsar ac¬ 
cording to the standard inner product definition: 

(S/N)^ = {h{X)\h{\)) = {h{XfG){G'^C{e)G)~^{G'^h{\)). 

(28) 

For each injected CGW, we draw the noise parameters of each 
pulsar {9a) from the corresponding posterior distribution of the 
SPA. We injected 1000 detectable CGWs, and each signal con¬ 
tained both Earth and pulsar terms. In figure[T]we plot the average 
(over the 1000 simulations performed) build-up of p = (S/N)^ as 
we add pulsars to the array, from the best to the worst. The plot 
shows the cumulative S/N square over the total S/N square, i.e., 
'Ylii-caP^-lP- ranking reveals that the array is heavily domi¬ 

nated by PSR J1909—3744, contributing more than 60% of p, fol- 
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lowed by PSR J1713+0747 at about 20%. We decided to form 
a dataset considering only the best pulsars building up 90% of 
p. We ended-up with 6 pulsars (PSRs J1909-3744, J1713+0747, 
J1744-1134, J0613-0200, J1600-3053, and J1012+5307) which 
constitute the restricted dataset mentioned above, and is the same 
used in our companion isotropic and anisotropic GW background 
searches ([Lentati et al.|2015||Taylor et al.|2015|l. 


4 RESULTS 

We now turn to the description of the main results of our analysis. 
As in the previous section, we present the outcome for the frequen- 
tist and Bayesian analyses separately, discussing first the detection 
results and then the upper limit computation for the various tech¬ 
niques. No evidence for CGW signals was found in the data; a sum¬ 
mary of all the 95% upper limits on the GW strain amplitude A as 
a function of frequency is collected in figure]^ 


4.1 Frequentist analysis 

4.1.1 Tp-statistic 

To tackle the issue of detection, we have evaluated IFp Nj = 99 
independent frequencies in the range [l/T, 2 x 10~^Hz] in steps 
of A/ = 1 /T, where T is the maximum observation time. In this 
dataset l/T = 2.0 x 10“®Hz. In the absence of a signal, 21Fp 
follows the PDF with n — 2N degrees of freedom given by 
equation l |18| >. The false alarm probability (FAP) associated to a 
given threshold Fq is simply given by the integral of the PDF and 
takes the form: 

roo n/2 —1 ^ 

PiFp>Fo)= / po{Fp)dFp = eM-J^o) V (29) 

To assess the global probability of finding a given value of Fp in 
our dataset, we need to take into account that we are evaluating it at 
99 different frequency bins, i.e., we are performing 99 independent 
trials. The global FAP is therefore given by 

FAP = 1 - [1 - P{Fp > Fo)ff . (30) 

We choose thresholds in Fq that correspond to FAPs of 0.01 and 
0.001. The results are presented in figure For the ML noise pa¬ 
rameters (solid line with circles and corresponding histogram on 
the right), the data is consistent with the noise description with a 
p-value of 0.93 and that there is no excess at any frequency above 
the 0.01 FAP threshold. However, the choice of the noise parame¬ 
ters, and hence the covariance matrix C, is crucial in evaluating Fp, 
and the SPA reveals that many parameters are poorly constrained. 
We therefore chose to create a whole distribution of Fp at each fre¬ 
quency, sampling over noise parameters from the posterior distri¬ 
butions produced by the SPA. This is overplotted as the yellow area 
in figure]^ The results obtained in this way are independent of the 
particular ML value and give a better view of the uncertainties in¬ 
volved. At each frequency, /, only pulsars with baselines T > 1// 
have been taken into account. This explains the rising FAP thresh¬ 
olds. The uncertainty in Fp induced by the uncertainty in the noise 
parameters is much larger at lower frequencies, where red noise and 
DM dominate the noise model. For the lowest frequencies, the ML 
evaluation of Fp does not even lie within the 90% region shown in 
the plot, which is a consequence of the fact that we sum the con¬ 
tributions of 41 broad, and often skewed, distributions. In several 




- - pdf 


— 1% FAP 



f [Hz] 


Figure 2. -statistics. The blue line represents Fp evaluated at 99 inde¬ 
pendent frequencies, for 41 pulsars, using the ML noise parameters. The 
right panel shows the histogram of the Fp values at all frequencies, and the 
dashed line is a central -distribution, which is the expected distribution of 
the Fp statistic in absence of a GW signal. The two are consistent with a p- 
value of 0.93, which is indicative of excellent agreement. The yellow area 
denotes the central 90% of Fp evaluated across the whole sample of noise 
parameters. The thresholds turn over below 6 nHz because of the reduced 
number of pulsars that have enough observing time span. 


pulsars the ML red noise or DM amplitudes tend to lie at the up¬ 
per end of the distribution, which leads to small values of Fp and 
results in this offset. 

After confirming that the SPA describes the data appropriately 
as only noise, we want to know how large a CGW contribution must 
be in order to make the -distribution non-central and clearly dis¬ 
tinguishable from the noise. This yields an upper limit on the GW 
strain A that might be present in our data and still consistent with 
the observed Fp value. A standard way to do this in frequentist 
analysis is through signal injections. We first fix the noise to the 
SPA-ML value; in this case the procedure to obtain the upper limit 
on A at each frequency / is as follows (see, e.g., [Ellis et al.|2012| l: 

(i) compute Fpp using the dataset; 

(ii) create 1000 different mock datasets i, by injecting in each of 
them one source with fixed strain A but otherwise random param¬ 
eters, and compute Fpp', 

(iii) compute the fraction y of mock datasets where Fpp > 

-Fp,o; 

(iv) repeat steps (ii) and (iii) with different A until y = 0.95. 

In practice, we iterate over a grid in logio.4 with step size 
0.1 and interpolate to find the point for which y = 0.95. We also 
want to obtain an upper limit that takes into account the uncertainty 
of the single-pulsar noise parameters. In doing this, the procedure 
outlined above is modified as follows: 

(i) compute a distribution of 1000 Fpp at each frequency / us¬ 
ing the dataset and 1000 noise parameters drawn from the posterior 
noise PDF obtained from the SPA; 

(ii) create 1000 different mock datasets i, by injecting in each of 
them one source with fixed strain A but otherwise random param¬ 
eters in each, and compute Fpp, each time drawing different noise 
parameters from the posterior noise PDF; 

(iii) compute the fraction y of mock datasets where Fpp > 
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Figure 3. 95% upper limit on the GW strain A obtained with the J^p- 
statistic. The blue line corresponds to noise parameters fixed to the ML 
values obtained in the time-domain SPA, while the dashed lines take into ac¬ 
count the full uncertainty in the noise estimation by sampling from the PDF 
distributions of either the pure time domain (green) or the time-frequency 
(red) SPA. At 1/year (the peak), the limit is poor because the GW signal at 
this frequency is absorbed in the fitting of the pulsar positions. 


tFp,o, where the single reference tFpfi is chosen to he the mean 
of the distribution obtained in step (i); 

(iv) repeat steps (ii) and (iii) with different A until y — 0.95. 

In other words, we want 95% of the tFp^i distribution to lie above 
the mean value tFpfi of the observed distribution. The motivation 
behind criterion (iii) is that it resembles the criterion for the ML 
upper limit and it is much more conservative than a Kolmogorov- 
Smirnov test. It is also possible to choose even more conservative 
criteria, susceptible to possible non-Gaussian tails in either distri¬ 
bution, but we do not explore this possibility here. 

Upper limits as a function of / are shown in figure]^ In all 
cases, the most stringent limit is around 6 — 7nHz, and reaches 
down to ^ = 6 X 10“^® for the Tp-ML case. We note, however, 
that when we allow the noise to vary, we get a limit which is a factor 
« 2 worse at low frequency, yielding a value of^ = 1.1 x 10“^^. 
This is consistent with figure at low frequency the ML estima¬ 
tor of the noise parameters is not representative of the noise pos¬ 
terior distribution, resulting in tFp values which are biased toward 
low values. Injections with lower A will therefore result in a de¬ 
tectable excess Fp, pushing down the upper limit. Note that our 
upper limit gets significantly worse at / < 3 x 10“®Hz, because 
red noise becomes significant for some pulsars, and not all 41 pul¬ 
sars in our array contribute down to those frequencies, having an 
observation time T < 10 years. At high frequency the 95% up¬ 
per limit degrades approximately linearly with /, consistent with a 
white-noise-dominated dataset. 


4.1.2 Fe-statislic 

Since Fe, is also a frequentist technique, the procedures to assess 
detection and to place upper limits are analogous to the Fp case. 
Here, in absence of signal, 2Fe, follows the PDF with n = 4 
degrees of freedom given by equation j23| >. The FAP associated 
with a given threshold To is simply given by the integral of the 



Figure 4. Result of the J^e-statistic detection analysis. The points are the 
trials of the search. The horizontal lines are the detection thresholds for 
different FAPs. 


PDF and takes the form: 


POO 

P{F, >Fo)= / po{Fe)dFe = (1 + Fo)e-^°. 
Jfo 


(31) 


Again, the global false alarm rate depends on the number of trials, 
according to equation l |30^ , which is now given by the number of 
independent templates in the sky location-frequency 3-D parameter 
space. The vast majority of templates we have used in the search are 
strongly correlated. We estimate the number of independent trials 
by constructing a stochastic 3-D template bank (see |Babak |2008t 
|Harry et al.|2009| l. We use a minimal match equal to 0.5 as the cri¬ 
terion of independence among different templates, and obtain 4276 
independent points in the searched parameter space (full sky and 
frequency band restricted to 2 — 400nHz). Figure presents the 
result of the detection analysis. The maximum of Fe, .Te,max, is 
found at f = 66 nHz and = {51.9°, 136.4°}. It corre¬ 

sponds to a FAP of 7%, which is compatible with a non-detection. 
Note that a GW signal with S/N= 5, if present in the data, would 
correspond to a FAP of « 5%, which is clearly too high to make 
any confident claim of detection, as shown in figure]^ 

Since the data are compatible with describing only noise, we 
can again compute the 95% upper limit on the strain amplitude A 
of a putative CGW as a function of / by means of signal injec¬ 
tions. The procedure is similar to the one employed in the Fp-ML 
analysis. Here, we construct a grid of frequencies, and at each grid 
point we consider a small frequency interval A/ = 1 nHz which is 
sufficient to capture the injected Earth term: 


(i) compute Fe,max,o, i-e. the maximum of Fe on the whole sky 
and in the narrow frequency band A/ in the raw dataset; 

(ii) create 1000 different mock datasets i, by injecting in each of 
them one source with fixed strain A but otherwise random param¬ 
eters, including both Earth and pulsar terms; 

(iii) for each dataset z, run a search (stochastic bank -l- multi- 
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Figure 5. The 95% upper limit on the GW strain obtained with the phase- 
marginalised Bayesian analysis, seai'ching on noise and signal simultane¬ 
ously {Bayes_EP_NoEv_noise, blue-circled curve), is compared to the same 
analysis with ML-fixed noise {Bayes_EP_NoEv, black curve) and to the 
noise-sampling J^p analysis (red curve). Details of each specific analysis 
are given in Table "TD" stands for time-domain analysis as opposed to 
the time-frequency approach developed in |Lentati et al.|j2013| . 

search genetic algorithnj^ to find i-e., the maximum of 

T'e on the whole sky and in the narrow frequency band A/; 

(iv) compute the fraction y of the mock datasets in which 

Te ,max,i > ^ 2,max,0. 

(v) repeat steps (ii) and (iii) increasing A until y — 0.95. 

The 95% sky-averaged upper limit obtained in this way is shown 
by the red curve in figure]^ and is in agreement with the results 
obtained with other methods. 


4.2 Bayesian analysis 

The frequentist analysis presented above already provided strong 
evidence against the presence of a signal in the data. Nonetheless, 
this can also be addressed in the Bayesian framework through the 
computation of the odds ratio defined by equation \21) . Since we 
give a priori no preference to the presence or absence of a signal 
in the data, we set the prior probability ratio to unity, and the odds 
ratio coincides with the Bayes factor. The Bayes factor is then sim¬ 
ply the ratio of the Bayesian evidence computed for the hypothesis 
Hi and hypothesis Ho as given by equation \21\ , which in our case 
reduces to: 

^ ^ j c{l\\5tMl\)ded\ 

/ c{e\5t)n{e)de 


A full use of the multi-modal genetic algorithm for each injection is 
computationally expensive and not needed, in practice. We therefore use 
a lighter and faster search for the injected signals. We construct a stochas¬ 
tic bank with minimal match 0.95 and we filter the data through this bank. 
We then identify the maximum of across the bank and refine our search 
running the genetic algorithm with 64 organisms evolved over 1000 gen¬ 
erations. The stochastic bank is generated only once for the full parame¬ 
ter space and contains 532 488 templates. In each search we use only the 
portion of template bank covering the parameter space region around the 
injected signal. 


In the case of fixed noise, we assume that the noise parameters 9 
are known exactly (fixed at their ML value), and the Bayes fac¬ 
tor is directly computed from the likelihood ratio multiplied by 
the priors, integrated over the source parameters A. We compute 
the evidence using both MultiNest and parallel tempering MCMC 
searches. In all the Bayesian searches with fixed noise we obtain 
Bayes factors close to zero, consistent with a non-detection and 
with the outcome of the frequentist analysis. In particular, we get 
log(I?) = —0.27 for the search Bayes_E, and log(S) = —0.31 for 
the search Bayes_EP_NoEv. 

The Bayesian analysis also returns samples from the joint 
posterior probability distribution of all model parameters. The 
marginalized distribution of any parameter of interest can then be 
evaluated by integrating (i.e., marginalizing) the joint posterior dis¬ 
tribution over all other parameters. We are particularly interested 
in the strain amplitude A. We can then split the vector parame¬ 
ter A = (.4, A') and integrate over A' to obtain the marginalized 
posterior for the parameter A. In practice, we divide the frequency 
range in small bins in which we carry out this marginalization pro¬ 
cedure separately. The 95% upper limit at each frequency corre¬ 
sponds to the value A for which 95% of the posterior distribution 
lies at .4 < .4; namely 

0.95= f dA 
Jo 

Results are shown in figure which compares all the upper 
limits on the GW strain achieved by all methods presented in this 
paper. For the non-evolving source case, the Bayes_EP_NoEv up¬ 
per limit agrees particularly well with the fixed-noise Fp-ML statis¬ 
tic. This is encouraging, since the two methods are similar in spirit 
as they adopt the same signal model and assume fixed/known noise 
parameters. For the evolving source case, the upper limit is 
very similar to both Tp-ML and Bayes_EP_NoEv, mimicking al¬ 
most perfectly their behavior at low frequency. The upper limits 
obtained by both the Bayes_E and the BayesJEP searches are nois¬ 
ier and slightly higher, but overall consistent with the others within 
a factor of two. 

As mentioned in sectionj^ we also ran a full 37-dimensional 
search over noise and signal parameters on the restricted set of the 
6 best pulsars in our PTA (c.f. figure [T](, and in a restricted fre¬ 
quency range of 5 — 15nHz where we have the best sensitivity. We 
used the phase-marginalised Bayesian analysis for non-evolving 
sources, and labeled the run Bayes_EP_NoEv_noise. The 95% up¬ 
per limit obtained in this case is shown in figure]^ together with 
the fixed noise Bayes_EP_NoEv and the noise-sampling Tp results. 
The Bayes_EP_NoEv_noise limit lies a factor 1.1 — 1.5 above the 
Bayes_EP_NoEv one. This is in line with the findings of |Arzouma-| 
|nian et aL] ( |2014^ , and confirms that our ML fixed noise upper limits 
are reliable within a factor < 1.5. It is also interesting to see that 
the Bayes_EP_NoEv_noise limit agrees fairly well with the noise¬ 
sampling Tp one. By analysing figure|^we can conclude that all the 
upper limits yielded by the different techniques agree within a fac¬ 
tor of two. We also observe that methods based on fixed noise (ML) 
parameters slightly underestimate the upper limit, which could be 
because the ML values are not always representative of the poste¬ 
rior distribution of the noise parameters. 


4.3 Sky maps 

Most of the searches outlined above are also sensitive to the source 
location on the sky. We can therefore extend our study and produce 


d\'C{A,\'\5t)n{A)ii{\'). (33) 
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Figure 6. The 95% upper limit on the gravitational wave strain for the 3 frequentist methods, i.e. tFp varying noise (Fp), Tp fixed noise (Fp_ML) and and 
the 3 bayesian methods, i.e. “evolving source” with Earth term only (Bayes_E) and with Earth and Pulsar terms (Bayes_EP) and “non-evolving source” with 
Earth and Pulsar terms (Bayes_EP_NoEv), see Table[T]for details. 


sky maps of the 95% upper limits provided by our analysis as a 
function of frequency. 

In the frequentist framework, this is straightforward to do in 
the context of the J^e-statistic, since it is sensitive to sky position 
(unlike tFp). As already mentioned, the upper limit is evaluated 
through massive signal injections according to the procedure out¬ 
lined in section l4.1.2l The difference is that now we have to divide 
the sky into “cells” and inject 500 sources at each cell location. This 
is much more computationally expensive than the evaluation of the 
sky-averaged limit; we therefore generate the sky map at 6.3nHz 
only, corresponding to our best sky-averaged limit. This is shown 
in figure As expected we are more sensitive in the region of the 
best pulsars. 

We can also produce targeted upper-limits as a function of 
sky location by means of Bayesian techniques. The problem here 
is that by splitting the posterior samples on a 3-D frequency-sky 
location grid, we end up with only a handful of points per cell, 
which are not enough to derive a reliable 95% upper limit. To 
mitigate this issue, we divided the sky into 4x2 patches on the 
(j> and 6 coordinates, respectively. A dedicated Bayesian analysis 
(fixed noise with marginalization of pulsar phases) on each patch 
yielded enough samples to sub-divide the region into a further 4x4 
sectors, for a total of 16x8= 128 resolution elements across the 
whole sky. Figure illustrates the sky map obtained in this way 
at a GW frequency of 7 nHz (a movie showing the evolution of 


the sky map across the relevant frequency range is available at 
http: // w w w.epta. eu. org/aom.html I. 

The qualitative agreement between the two maps is quite 
good. In both of them, the best pulsars are shown as white dots, 
with size proportional to their contribution to the square of the S/N. 
As expected, the most constraining (i.e. lowest) upper-limits on the 
strain of a putative CGW lie around the location of the best pul¬ 
sars in the array, and the sky maps shows a clear dipolar pattern. 
The closest galaxy clusters in the Universe, i.e. Virgo and Coma, 
are located at the transition between the two regions of the dipolar 
pattern, in an area of “average sensitivity”. 


5 ASTROPHYSICAL INTERPRETATION 

The upper limits on CGWs from SMBHBs presented in this paper 
are currently the most stringent in the literature. We turn now to 
investigate their impact on the astrophysics of SMBHBs. 


5.1 Horizon distance 

Each of the 95% upper limits on A derived in the previous section 
can be easily converted into a horizon distance for CGW detection 
as a function of mass and frequency using equation 
is the strain upper limit as a function of frequency obtained with a 
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Figure 7. Sensitivity sky map at / = 6.3 nHz computed with the Te. com¬ 
puted with 500 injections in 48 directions in the sky ("cells"). The color 
scale corresponds to log^^Q of the 95% upper limit on the strain amplitude 
A. The white points indicate the positions of the 6 best pulsars with sizes 
corresponding to their contribution to the S/N. Black dots indicate the loca¬ 
tion of the Virgo and the Coma clusters. 



Figure 8. Sensitivity sky map at / = 7 nHz computed with the phase- 
marginalised Bayesian technique for a “non-evolving source”. The color 
scale and points are the same as for figure]?] 


specific method, then 

a^5/3 

= (34) 

In a frequentist sense, this has to be interpreted as the distance at 
which, on average, a source of mass Me emitting at frequency / 
located anywhere on the sky would result in a value of the detection 
statistics higher than what we measure in the data with 95% proba¬ 
bility, if it was there. As an example, results for the Tp-ML statistic 
are presented in figure An interesting feature of the plot is that, 
for a given Me, Dh is essentially constant (slowly declining) for 
/ > 5 X 10“®Hz. This is because of the cancellation effect between 
the rising CGW amplitude with frequency, A oc and the 

PTA sensitivity, which degrades almost linearly with / (see figures 
[^and|^. In this frequency range, and with the current sensitivity, 
we can exclude the presence of a SMBHB with Me > IO^Mq 
out to a distance of about 25Mpc, i.e. well beyond the distance to 
the Virgo cluster, and with Me > 3 x 10® Mq out to a distance 
of about 200Mpc, i.e. twice the distance to the Coma cluster. Note 
that Virgo and Coma themselves are located in a region of “average 
sensitivity” in our sky sensitivity map (see figure]^, meaning that 
we can rule out the presence of SMBHBs (with the characteristics 
described above) in these specific clusters. We remind the reader 
that these numbers are for SMBHBs with a given redshifted chirp 



Figure 9. Horizon distance as a function of GW frequency for selected 
values of Me, based on the Tp-ML upper limit. 


mass; because of the 1 + z factor, horizon distances evaluated at the 
same values of the intrinsic SMBHB mass will be slightly larger. 

A number of potentially interesting sources have been pro¬ 
posed in the literature. Among these are the two blazars OJ287 
UValtonen et aL]|2008t and PG 1302-102 ( [Graham et aT][20T5l ). 
Both objects are located at z ~ 0.3 corresponding to Dl ~ 
l.SGpc. At such distances, we cannot rule out any system below 
Me ~ IO'^^Mq, which far exceeds the plausible range of chirp 
masses inferred for these two objects. Should their binary nature be 
confirmed, these two systems would likely be among the thousands 
of contributors to the stochastic GW background. At this stage, 
SMBHBs need to be more massive and/or nearby to be resolved 
by a PTA. 


5.2 Probability of detection 

A natural question that arises at this point is: could we expect a 
detection of a CGW signal with the current EPTA sensitivity? We 
now evaluate the probability of detecting CGWs from an individual 
SMBHB with an array like the current EPTA, by using a large set 
of observationally-based simulations of the cosmic SMBHB pop¬ 
ulation. Each simulation represents a particular realisation of the 
ensemble of SMBHBs. In a nutshell, galaxy merger rates are ob¬ 
tained using a selection of galaxy mass functions and close galaxy 
pair fractions from the literature; merging galaxies are populated 
with SMBHs following empirical black hole-galaxy host relations; 
finally, each binary is assumed to emit GWs while inspiralling in 
a quasi-circular orbit. Given the broad range of different models 
taken into account and their uncertainties, numerous simulations 
are created in order to cover all possible configurations consistent 
with the observations. The SMBHB populations obtained in this 
way are consistent with the results of semi-analytic halo merger 
trees and cosmological N-body and hydrodynamical simulations. 
More details on the simulations and the models employed to pro¬ 
duce them can be found in |Sesan^j2013| l. 

One can perform signal injections drawing the sources from 
these models and run all the different detection pipelines detailed 
above, to assess detection probabilities. However, this is an expen¬ 
sive task, and we do not need such a refined analysis at this stage. 
We instead simplify the problem following a similar approach as in 
[Rosado et al.| ( [20T5| l. For a given realization of the SMBHB pop¬ 
ulation, we group GW sources in frequency bins A/ = l/T, and 
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Figure 10. GW strain amplitude versus GW observed frequency. The 
coloured lines represent the different upper limits presented in this work. 
The shading gives the probability of detecting a SMBHB in a particular in¬ 
terval of strain and frequency. That detection probability increases towards 
lower frequencies and smaller values of strain (on the lower-left corner). In 
the legend, the percentage of detection probability is given for each of the 
upper limits. 


compute the characteristic strain 


h 


2 

c 


J2k ^kfk 

A/ 


(35) 


the sum runs over all binaries falling in the frequency bin. We then 
identify the loudest source in each bin, and compute its S/N follow¬ 
ing [Sesan^^^^cchl^ ([20T^ assuming that the noise is given by 
the sum of the strains of the GWs produced by all other binaries. 
In practice, we are assuming that all other sources produce an “un¬ 
resolved background”, and we check whether the loudest source 
“sticks out” of it. We assume a detection statistic described by a 
distribution with 4 degrees of freedom, and we consider “individ¬ 
ually resolvable” only those sources with S/N surpassing the FAP 
threshold of 0.1% related to this distributior0 

Let us assume a particular upper limit on the GW strain am¬ 
plitude, among those presented in figure]^ and call it UL®. At a 
particular frequency bin fj , we simply estimate the probability of 
detecting a SMBHB with such sensitivity as the fraction of real¬ 
isations in which a resolvable binary produces a strain amplitude 
A > Ag 5 %{fj). We call this detection probability p(D|UL®, fj). 
The probability of not detecting a binary at that frequency is thus 
p(N|UL®, fj) = 1 — p(D|UL®, fj). Assuming that the probabili¬ 
ties of different frequency bins are independent, the probability of 
not detecting a binary in any frequency bin is the product of the 
individual values p(N|UL®, fj). Hence, 


p(D|UL®) = 1 - p[(l - p(D|UL®, fj)) (36) 

3 


is the probability of detecting a SMBHB at any frequency bin, for 
the upper limit UL®. 

The detection probability at any frequency obtained for each 
of the upper limits is given in the legend of figure[^ The maximum 
detection probability achieved with the EPTA upper limits is below 
~ 1%. Therefore, we can safely conclude that a non-detection is 
consistent with the theoretical expectations. 


® By assuming a distribution with 4 degrees of freedom, we are assum¬ 
ing Te as detection statistic. This is an arbitrary choice dictated by compu¬ 
tational convenience only. Results are, however, qualitatively unchanged if 
a different statistic (e.g. Fp) is assumed. 



Number of frequency bins shifted 



Number of frequency bins shifted 


Figure 11. Cumulative PDF of the frequency shift between pulsar and Earth 
terms, in units of the adopted frequency bin A/ = 15yr~^. The upper and 
the lower panels assume a pulsar distance of 2kpc (average distance of the 
EPTA pulsars) and 6kpc (maximum distance), respectively. 


5.3 Frequencies of the Earth and the pulsar terms 

In our searches, we distinguished between evolving and non¬ 
evolving GW signals, presenting distinct search methods for each 
of them. One may therefore ask, whether one type of signal is more 
likely than the other, in order to better focus development efforts on 
specific analysis pipelines. We can use the same simulated SMBHB 
populations discussed above to answer this question. As shown by 
figure [To] only a small percentage of them leads to a detection with 
the sensitivity of the current EPTA, being therefore inconsistent 
with observations. Nevertheless, it is still meaningful to study the 
outcome of those realizations, since they would resemble the true 
ensemble of SMBHBs in the fortunate case of a detection in the 
near future. 

For each of the observable SMBHBs in those realizations, we 
calculate, according to equation ([^, the frequency evolution of the 
emitted GW after a time lapse of 4kpc/c, which would be the max¬ 
imum time difference between pulsar and Earth terms for a pulsar 
located at 2kpc (which is approximately the mean distance to the 
EPTA pulsars). Erequency shifts are defined as the difference be¬ 
tween the GW frequency before and after that time lapse in units 
of the frequency resolution bin of the array (assumed to be 15yr“^ 
here). Their distribution is shown on the upper plot of figure m 
The lower plot is analogous, but assuming a pulsar located at 6kpc 
(which corresponds to approximately the largest distance to a pul¬ 
sar in the EPTA). Should EPTA detect an individual source in the 
near future, it could either be evolving or non-evolving with nearly 
equal probability. Even considering the shifts produced in the fur¬ 
thest pulsar only (5/ 6kpc), there is still a 36% probability that 
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the source would be non-evolving. Decreasing the PTA sensitiv¬ 
ity floor would make it sensitive to lower mass binaries (which 
evolve faster), but would also improve the chance of detection at 
higher frequency, where evolution is more likely. Likewise, extend¬ 
ing the observation time will allow to see binaries at lower fre¬ 
quency, where evolution is less likely; it will, however, also shrink 
the size of the frequency bin (A/ = l/T), making it easier for 
a source to sweep through different resolution elements. Detection 
strategy development for both classes of sources is therefore war- 
ranteqfl 


6 CONCLUSIONS 

In this paper, we searched for continuous gravitational wave signals 
in the latest release of the European Pulsar Timing Array dataset. 
We adopted both frequentist and Bayesian techniques, searching 
for both frequency-evolving and strictly monochromatic signals. In 
most of the cases, we fixed the value of the noise parameters in each 
pulsar to the maximum likelihood estimated in a separate single¬ 
pulsar analysis. This choice was primarily dictated by computa¬ 
tional feasibility, but is certainly non optimal, since pulsar noise 
and GW signals might be degenerate and one should include both 
simultaneously in the search. To validate our results we therefore 
also performed a frequentist analysis, sampling from the posterior 
distributions of the noise parameters returned by the SPA (simply 
labeled and a full Bayesian search over both the noise and 
the signal (labeled Bayes_EP_NoEv_noise). Because of the high 
dimensionality of the search parameter space, the latter has been 
conducted on a restricted dataset including only the 6 best pulsars 
in the array. 

None of the analysis yielded any evidence of the presence of a 
signal, and only upper limits on the amplitude ^ of a putative CGW 
could be placed. The excellent quality and length of the dataset al¬ 
lowed us to place limits comparable to those of |Zhu et al.| j2014^ 
at / > lOnHz and a factor of two better at / < lOnHz, yeld- 
ing the overall most stringent constrains to date. All the employed 
methods yield 95% upper limits on A (. 495 %) consistent within 
a factor of two across the whole 2nHz-400nHz frequency range. 
Our best sensitivity is in the 5nHz-7nHz interval, where we find 
6 X 10“^® < . 495 % < 1.5 X 10“^^, depending on the adopted 
method. The most robust analysis {Bayes_EP_NoEv_noise) results 
in . 495 % = 9 X 10“^® at 6 nHz. Limits on the strain amplitude 
can be converted to horizon distances as a function of source mass 
and frequency. We exclude the existence of SMBHBs with sepa¬ 
ration < O.Olpc and Afc > 10 ®Mq out to a distance of about 
25Mpc (well beyond Virgo), and with Me > 1 O® ®M 0 out to a 
distance of about 200Mpc (twice the distance to Coma). In recent 
years, several “overmassive” black holes have been found in the 
local Universe, with measured masses in excess of 10 ^®Mq. Our 
analysis excludes that any such system lives in a compact binary 
within a distance of about IGpc {z ~ 0.2). Finally, we compared 
our limits to the predictions of state of the art models of the cosmic 
population of SMBHBs. We found a detection probability of < 1% 
at current sensitivity, consistent with the null result of our searches. 

The present analysis has also highlighted a few interesting 
technical issues related to the search methods and to the nature 
of the dataset. Despite not being robust for detection purposes, as 

® A more detailed study of the expected properties of the first detectable 
SMBHBs can be found in|Rosado et al.|j2015^. 


pointed out by | Arzoumanian et al.| ( |2014) , fixed noise analysis up¬ 
per limits are consistent within 50% of those obtained by searches 
over the full parameter space (i.e., including signal and noise si¬ 
multaneously). Therefore, so long as the data do not support the 
presence of a signal, a computationally cheap analysis of this type 
can be carried out over an extensive dataset of numerous pulsars, 
possibly yielding more interesting astrophysical constraints on the 
low-redshift SMBHB population in the near future. Eventually, si¬ 
multaneous searches over the signal and noise parameters will be 
required for a confident detection claim. However, those are ex¬ 
tremely expensive, and novel techniques capable of efficiently han¬ 
dling parameter spaces of lOO-l- dimensions must be developed. The 
reason why the results of the full search on the restricted dataset of 
6 pulsars is consistent with those provided by fixed noise analy¬ 
sis on the full set of 41 pulsars, is that the current EPTA array is 
heavily dominated by a handful of ultra-stable MSPs. In particular, 
PSRs J1909—3744 and J1713-I-0747 combined account for 80% of 
the EPTA sensitivity to CGWs. As a result, the EPTA dataset sen¬ 
sitivity has a strongly dipolar pattern across the sky, varying by 
almost a factor of four over the celestial sphere. The discovery of 
new ultra-stable MSPs will therefore be crucial to provide a better 
sky coverage, ensuring that no “blind spots” are left, and thus en¬ 
hancing the probability of detecting CGWs in the coming decade. 
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